Introduction

Nous sommes deux étudiants habitués aux métros parisiens, surchargés tout autant que le trafic routier, comme dans d’autres grandes villes. Nous pouvons donc affimer que se déplacer en vélo présente beaucoup d’avantages…mais malheureusement, un inconvénient existe : parfois il n’y a plus de vélos disponibles dans une station désirée, ou alors plus de places pour les rendre.

Nous avons lu un article d’un enseignant chercheur en mécanique des fluides qui a créé justement un logiciel de prédiction, intégré dans l’application «La Bonne Station» de Keolis, destinée aux utilisateurs du VCub a Bordeaux (équivalent des velibs a Paris) et nous avons trouvé intéressant cet enjeux de prévision. Son objectif est de maximiser l’efficacité des opérations de régulation des vélos: Quelles stations approvisionner en priorité ? Quelles stations agrandir en priorité ? Où positionner de nouvelles stations ? Quelles stratégies adopter pour être le plus efficace possible ?

Présentation des données

Pour notre projet d’analyse prédictive reliée au vélibs, nous allons travailler sur le jeu de données London bike.
Nous l’avons trouvé sur le site Kaggle. Il correspond a des données de locations de vélos dans la ville de Londres.

Appercu du jeu de données
timestamp cnt t1 t2 hum wind weather holiday weekend season
2015-01-04 00:00:00 182 3.0 2.0 93.0 6.0 3 0 1 3
2015-01-04 01:00:00 138 3.0 2.5 93.0 5.0 1 0 1 3
2015-01-04 02:00:00 134 2.5 2.5 96.5 0.0 1 0 1 3
2015-01-04 03:00:00 72 2.0 2.0 100.0 0.0 1 0 1 3
2015-01-04 04:00:00 47 2.0 0.0 93.0 6.5 1 0 1 3

Il contient \(17414\) observations, recensées sur deux ans, pendant chaque heure entre le 04/01/2015 à 00h00 et le 03/01/2017 à 23h00, et \(10\) variables :

\(\underline{Problématique}\)
Nous commencerons par chercher a prédire le nombre d’emprunts de vélos des 4 derniers mois du jeu de données, avec des modèles qui utilisent juste les variables explicatives.
Puis nous ferons un modele prédictif, afin de prevoir les futurs emprunts de vélos.

Notre variable d’intérêt sera donc le nombre de vélos empruntés «cnt».
Et nos variables explicatives seront les autres.

\(\underline{Traitement~des~données}\)
Nous commençons par verifier que nous n’avons pas de données manquantes:

C’est bien le cas, donc nous n’avons pas besoin de faire d’imputation de celles-ci!

Ensuite, nous ajoutons plusieurs variables qui nous seront utiles pour la suite:

Appercu du nouveau jeu de données
timestamp cnt t1 t2 hum wind weather holiday weekend season day_of_month month day_of_week hour Posy date Year
2015-01-04 00:00:00 182 3.0 2.0 93.0 6.0 3 0 1 3 4 1 1 0 0.0082418 2015-01-04 2015
2015-01-04 01:00:00 138 3.0 2.5 93.0 5.0 1 0 1 3 4 1 1 1 0.0082418 2015-01-04 2015
2015-01-04 02:00:00 134 2.5 2.5 96.5 0.0 1 0 1 3 4 1 1 2 0.0082418 2015-01-04 2015
2015-01-04 03:00:00 72 2.0 2.0 100.0 0.0 1 0 1 3 4 1 1 3 0.0082418 2015-01-04 2015
2015-01-04 04:00:00 47 2.0 0.0 93.0 6.5 1 0 1 3 4 1 1 4 0.0082418 2015-01-04 2015

Analyse descriptive

Voici le nombre de vélos empruntés au cours du temps:

Nos données ne recouvrent que deux ans, mais nous pouvons voir notamment une saisonnalité.
Nous n’avons pas de changement signification d’une année à l’autre dans la locations des vélos.
On remarque a priori qu’ils sont plus utilisés pendant l’été que l’hiver.
Nous voyons deux valeurs anormales / deux pics élévés durant l’été 2015. Les deux pics correspondent à des jours de grève des transports.

Regardons tout d’abord le nombre de vélos loués en fonction des différentes variables.

\(\underline{Variables~calendaires}\)

\(\bullet~Graphiques~annuel\)

\(\bullet~Graphiques~hebdomadaire\)

\(\bullet~Graphiques~quotidien\)

\(\underline{Variables~liées~a~la~météo}\)

Nous remarquons que les vélos sont le plus loués lorsqu’il fait plus chaud, et qu’il y a moins d’humidité.
La vitesse du vent n’a pas vraiment l’air d’influer sur les locations par contre.

\(\underline{Corrélations}\)

\(\underline{Autocorrélations}\)

Si la série était stationnaire, nous verrions pratiquement toutes les lignes dans les intervalles de confiance a 95% représentés en verts. Or, nous voyons que chaque pic est hors de ces lignes.

Modélisation

Predictions de la fin du jeu de données

Notre objectif est de développer des modèles basés sur l’ensemble de données d’entraînement et de faire des prédictions du nombre horaire de location de vélos à l’aide de l’ensemble de données de test.
Nous allons donc estimer des modèles sur des données d’entrainement et faire des prévisions sur des données de test.
Nous commençons par divisé notre jeu de données en une partie pour entrainer nos modèles, et une partie pour les valider.
Nous gardons les 4 derniers mois pour nos validations.

Nous avons 14488 observations pour nos tests, et 2926 pour nos validations.

Nous allons appliqué plusieurs modèles afin de prédire le nombre total de locations de vélos par heure, pour les quatres dernier mois.
Nous les évaluerons grâce a differentes valeurs :

  • La racine de l’erreur quadratique moyenne (RMSE) définie par: \(RMSE = \sqrt{MSE} = \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i–\hat{y}_i)^2}\)
    Elle s’exprime dans la même unité que la variable à prédire.
    Plus elle est élevée, moins le modèle est performant.
  • L’erreur absolue moyenne en pourcentage (MAPE) définie par: \(MAPE\) = \(\frac{1}{n} \sum_{i=1}^n \displaystyle\left\lvert \frac{e_i}{y_i} \right\rvert\) = \(\frac{1}{n} \sum_{i=1}^n \displaystyle\left\lvert \frac{ (y_i–\hat{y}_i) }{ y_i } \right\rvert\)
    Elle s’exprime en pourcentage et elle peut se traduire par le pourcentage moyen d’écart entre la valeur prédite et la réalité.
    De même, plus elle est élevée, moins le modèle est performant.
  • Le temps d’exécution

Remarque:
- Nous avons choisi d’utiliser la MAPE et non la MAE (erreur absolue moyenne), pour deux raisons:
La première est qu’un pourcentage est plus facile à comprendre, et la seconde est que les scores MAPE peuvent être comparés entre différents modèles. Cependant, si vos valeurs réelles auraient été proches de 0, nous aurions choisis la MAE!!
- La MSLE (erreur logarithmique quadratique moyenne) réduit l’importance des erreurs pour de grandes valeurs réelles, comme le fait deja la MAPE. Elle ne s’exprime pas dans une unité simple, ce qui la rend difficile à interpréter, ce pourquoi nous ne la choisissons pas non plus!

\(\underline{GLM}\)

RMSE = 1005.917 
MAPE = 1.536 
utilisateur     système      écoulé 
          0           0           0 

Ce modèle est très rapide, mais le résultat n’est pas terrible..

\(\underline{SARIMA}\)

Series: bikeCount 
ARIMA(5,1,0) 

Coefficients:
         ar1      ar2      ar3      ar4      ar5
      0.3228  -0.3496  -0.0672  -0.0203  -0.0606
s.e.  0.0076   0.0080   0.0084   0.0079   0.0076

sigma^2 = 416609:  log likelihood = -137366.9
AIC=274745.8   AICc=274745.8   BIC=274792.4

RMSE = 1328.121 
MAPE = 4.406 
utilisateur     système      écoulé 
       0.01        0.00        0.01 

Le modèle SARIMA est une extension du modèle ARIMA.
Nous otenons une RMSE d’environ 1328 vélos, ce qui est encore moins bon que pour notre modele de regression. On constate que le modèle est très bon sur les données d’entraînement et moins bon sur les données de test. On peut donc soupçonner des cas de sur-apprentissage.

\(\underline{Arbres~de~décision}\)

RMSE = 534.295 
MAPE = 0.788 
utilisateur     système      écoulé 
       0.24        0.02        0.25 
  • Sur le premier graphique, nous pouvons visualiser l’arbre obtenu.
    Trois de nos variables hour, t1 et weekend ont été utilisées, et 14 nœuds terminaux ont été générés.
    Sur notre ensemble d’entrainement, nous voyons que la variable d’heure hour est la première variable choisie pour la décision. Le nombre de vélos empruntés attendu, lorsque hour<7 est de 197, et 29% des échantillons tombent dans cette feuille. Lorsque hour est inférieure à 7heure, moins de vélos sont utilisés, et lorsque hour est >7h et <20h, la température commence à entrer en jeu.

  • Sur le deuxieme graphique, nous pouvons voir l’erreur de validation croisée (relative) pour chaque sous-arbre, du plus petit au plus grand.
    Dans “les coulisses”, rpart() applique automatiquement une validation croisée pour élaguer l’arbre. L’axe x en bas correspond a la complexité de l’arbre.
    L’axe x en haut correspond au nombre de nœuds terminaux, donc la taille de l’arbre.
    L’axe y a l’erreur de convergence.
    Un bon choix de cp pour l’élagage est souvent la valeur la plus à gauche pour laquelle la moyenne se trouve sous la ligne horizontale. Nous voyons que l’endroit optimal pour tailler l’arbre correspond a 12 nœuds.

  • Notre RMSE pour l’ensemble de test est d’environ 534 vélos!!
    la MAPE vaut environ 0.79, et le temps d’execution est de 0.118 secondes.

  • Remarque: Les arbres de décision sont des algorithmes où les hyperparamètres sont essentiels en ce qui concerne leur performance. Par exemple, un nombre de profondeurs maximales (maxdepth) trop élevées peut entraîner un surajustement ou une variance élevée, et si ce nombre est trop faible cela peut entraîner a l’inverse un manque d’ajustement ou un biais élevé. Il faut donc faire attention car l’objectif est d’avoir les meilleurs resultats, mais l’arbre pourrait être suradapté.

\(\underline{GAM}\)

RMSE = 603.294 
MAPE = 1.056 
utilisateur     système      écoulé 
       0.04        0.00        0.05 

Plus rapide (seulement 0.014 secondes), mais la RMSE est d’environ 603 vélos, et la MAPE d’environ 1.06 : c’est un peu moins bien que l’arbre !

\(\underline{Randomn~forest}\)

RMSE = 557.427 
MAPE = 0.835 
utilisateur     système      écoulé 
       0.16        0.11        0.14 

Notre RMSE pour l’ensemble de test est d’environ 549 vélos, ce qui est un tout petit peu plus que pour l’abre de décision, mais ils ont la même valeur de MAPE, environ 0.79, et nos forets aleatoires sont plus rapide (0.028 secondes)

\(\underline{XGBoost}\)

RMSE = 372.815 
MAPE = 0.718 
utilisateur     système      écoulé 
       0.04        0.00        0.02 

La RMSE vaut environ 373 vélos et la MAPE environ 0.718!!! Nous obtenons les meilleurs resultats, l’extreme gradient boosting dépasse tout le monde! De plus le temps d’execution est seulement de 0.002 secondes.

\(\underline{Perceptron}(MLP)\)

Entrainons le perceptron multi-couches en faisant varient le nombre de neurones, le nombre de couches afin de trouver un modèle qui prédit mieux.

La fonction RELU \(f(x)=x^{+}=max(0,x)\) sera notre fonction d’activation. On utilisera le package Keras.

  • cas d’une couche en entrée et sortie (pas de couche cachée).

128 neurones à la première couche et 1 à la dernière.

RMSE = 909.46 
MAPE = 1.345 
utilisateur     système      écoulé 
       0.44        0.82        0.32 
  • cas d’une couche en entrée, une couche cachée et une couche de sortie.

500 neurones à la première et deuxième couche et 1 à la dernière.

RMSE = 669.07 
MAPE = 0.542 
utilisateur     système      écoulé 
       0.67        1.03        0.53 
  • cas d’une couche en entrée, deux couches cachées et une couche de sortie.

600 neurones à la première,deuxième et troisième couche et un à la dernière.

RMSE = 635.375 
MAPE = 0.364 
utilisateur     système      écoulé 
       1.19        1.05        0.59 

On constate que le meilleur modèle est le 3 en termes de RMSE et MAPE. Donc plus le modèle aura des couches et de neurones mieux sera les résultats.

Prediction du futur après le jeu de données

Nous allons maintenant faire un modele prédictif pour essayer de prevoir le jour même, les 24 points (1 point par heure) pour le lendemain, avec nos données allant jusqu’à la veille.

\(\underline{Prophet}\)

Prophet se base sur un modèle additif où les tendances non linéaires correspondent aux saisons annuelles, hebdomadaires et quotidiennes, ainsi qu’aux effets des vacances.
Le modèle fonctionne mieux avec les séries temporelles qui ont de forts effets saisonniers et les données historiques de plusieurs saisons.
Il examine les données antérieures afin de pouvoir analyser le futur.
Il utilise un modèle de série temporelle qui peut être décrit par 3 composantes principales du modèle: les tendances, les saisons et les vacances, selon l’équation suivante :

\(y(t) = g(t) + s(t) + h(t) + \epsilon(t)\)

\(g(t)\) : changements non périodiques des valeurs de la série \(s(t)\) : changement périodique (saisonnier, hebdomadaire, annuel) \(h(t)\) : l’effet des jours fériés sur le calendrier, pouvant créer des irrégularités \(\epsilon(t)\) : changements inhabituels non pris en compte par le modèle

Prophet utilise des régresseurs de température et d’humidité. Tout cela nous semble très interessant dans notre cas.
Nous allons donc faire des prédictions pour le prochain jour pour chaque heure (ou pour 24 heures).

                       ds      yhat
17414 2017-01-03 23:00:00  158.7483
17415 2017-01-04 00:00:00  766.0758
17416 2017-01-04 01:00:00 1293.9043
17417 2017-01-04 02:00:00 1499.2879
17418 2017-01-04 03:00:00  301.4787
17419 2017-01-04 04:00:00 -115.5142
17420 2017-01-04 05:00:00  798.4119
17421 2017-01-04 06:00:00 1346.4293
17422 2017-01-04 07:00:00 1761.7837
17423 2017-01-04 08:00:00 3551.3491
17424 2017-01-04 09:00:00 2239.1037
17425 2017-01-04 10:00:00 1764.9763
17426 2017-01-04 11:00:00 1911.5235
17427 2017-01-04 12:00:00 1076.4066
17428 2017-01-04 13:00:00 2219.9535
17429 2017-01-04 14:00:00 1773.2393
17430 2017-01-04 15:00:00 1525.1560
17431 2017-01-04 16:00:00 2210.3772
17432 2017-01-04 17:00:00 3270.4694
17433 2017-01-04 18:00:00 3459.5608
17434 2017-01-04 19:00:00 2913.8282
17435 2017-01-04 20:00:00 2161.0340
17436 2017-01-04 21:00:00 1736.3712
17437 2017-01-04 22:00:00 1255.1763
17438 2017-01-04 23:00:00  635.5908

Nous retrouvons la prédiction de nos 24 points pour le lendemain.

Sur ce graphique les points en noirs representent les données réelles, et les lignes bleues sont les valeurs prédites. Nous le dernier pic bleu correspond a la prediction future de la journée suivante.

Nous voyons les performances du modèle grâce au nuage de points, entre la valeur prédite et la valeur réelle.
La droite ajustée (en rouge) a une pente positive, ce qui signifie que plus les valeurs “actual” sont grandes, plus les valeur “pred” le sont aussi.
Adjusted R-squared: 0,6261
[R-squared est un nombre compris entre 0 et 1, et plus il est grand, meilleur sera le modèle résultant] Notre résultat n’est pas le plus optimal.

Nous terminons cette partie modelisation par un test ARMA avec Prophet.
Nous allons essayer de faire une rediction du futur sur un an avec la moyenne mobile.

Aggregation

À présent agrégions les différents modèles obtenu jusqu’à présent, à savoir l’arbre de décision, le modèle additif généralisé, la forêt aléatoire, le gradient boosting et le perceptron multi-couches. On utilisera le package opera pour l’agrégation..

[1] 2926    5
Aggregation rule: EWA 
Loss function:  squareloss 
Gradient trick:  FALSE 
Coefficients: 
     rpart gam        rf xgboost mlp
 6.35e-216   0 5.02e-253       1   0

Number of experts:  5
Number of observations:  2926
Dimension of the data:  1 

        rmse  mape
EWA      373 0.716
Uniform  433 0.518

On peut toute de suite constater que l’algorithme n’a utilisé que les meilleurs modèles. Le gradient boosting domine tous les autres experts. Le perceptron multi-couches apparaît tout au début en éclipsant les autres experts.

Le RMSE et le MAPE du modèle agrégé sont de 373 et 0.715. Soit le meilleur mape que ceux de tous nos modèles précédents et le deuxième meilleur rmse derrière xbgboost de rmse 372.

Conclusion

Nous remarquons que pour certains modèles, les prévisions du nombre de vélos prennent parfois des valeurs négatives, ce qui n’est pas possible mais comme ce sont des prévisions, cela fait partie de “l’erreur”

Resumé des modèles utilisés
RMSE MAPE Temps d’execution
GLM 1005.91704439951 1.53614602226818 0.07
SARIMA 1328.12067903949 4.40586219999758 0.02
Arbre de décision 534.294666099615 0.787693158684004 0.118
Modèle Additif Généralisé 603.293713663435 1.05554785178295 0.007
Foret Aléatoire 557.427418914906 0.835415619899165 0.029
eXtreme Gradient Boosting 372.81500435415 0.717676653923912 0.002
MLP 635.374503996781 0.363700834869851 0.67
EWA 373 0.717 NA
Uniform 575 0.498 NA

Nous obtenons des MAPE < 10%, ce qui est très bien !

\(\underline{Validation~du~modèle~final}\)
Nous pouvons constater que c’est avec le modèle XGBoost que nous obtenu les meilleurs résultats, et c’est celle qui a également nécessité le moins de temps de calcul, car elle demande moins de paramètres à former que les autres modèles.

Les algorithmes XGBoost et Random Forest souffrent d’une variance élevée, alors on pourrait ajouter plus de données et augmenter l’effet de la régularisation… Nous avons entrainé nos modèles sur des données d’un an et 8 mois, et validé sur un ensemble de test de 4 mois. Les vélibs sont moins loués pendant l’hiver que pendant le reste de l’année. Peut-être que l’ajout de plus de données pendant l’hiver aidera à améliorer les performances pendant cette période, ou faire un sur-échantillonnage.

En conclusion générale nous sommes ravis d’avoir pu tester ces modèles sur un jeu concret de donnée.